Network Models

Modeling relationships and interactions between entities as a graph of nodes and edges.

A network represents the relationships (links) between entities (nodes). These links can be weighted (weighted network) or unweighted (binary network), directed (directed network) or undirected (undirected network). Regardless of their type, networks generate links shared by nodes, leading to data dependency when modeling the network. One proposed solution is to model network links with random intercepts and slopes. By adding such parameters to the model, we can account for the correlations between node link relationships.

Considerations

Caution

• The particularity here is that varying intercepts and slopes are generated for both nodal effects 🛈 and dyadic effects 🛈. These varying intercepts and slopes are identical to those described in previous chapters and will therefore not be detailed further. Only the random-centered version of the varying slopes will be described here.

Example

Below is an example code snippet demonstrating a Bayesian network model with a sender-receiver effect. This example is based on Ross, McElreath, and Redhead (2024).

Python

# Setup device------------------------------------------------
from BI import bi, jnp

# Setup device------------------------------------------------
m = bi(platform='cpu')
# Simulate data ------------------------------------------------
N = 50
individual_predictor = m.dist.normal(0,1, shape = (N,1), sample = True)

kinship = m.dist.bernoulli(0.3, shape = (N,N), sample = True)
kinship = kinship.at[jnp.diag_indices(N)].set(0)

def sim_network(kinship, individual_predictor):
  # Intercept
  alpha = m.dist.normal(0,1, sample = True)

  # SR
  sr = m.net.sender_receiver(individual_predictor, individual_predictor, s_mu = 0.4, r_mu = -0.4, sample = True)

  # D
  DR = m.net.dyadic_effect(kinship, d_sd=2.5, sample = True)

  return m.dist.bernoulli(logits = alpha + sr + DR, sample = True)


network = sim_network(m.net.mat_to_edgl(kinship), individual_predictor)

# Predictive model ------------------------------------------------

m.data_on_model = dict(
    network = network, 
    dyadic_predictors = m.net.mat_to_edgl(kinship),
    focal_individual_predictors = individual_predictor,
    target_individual_predictors = individual_predictor
)


def model(network, dyadic_predictors, focal_individual_predictors, target_individual_predictors):
    N_id = network.shape[0]

    # Block ---------------------------------------
    alpha = m.dist.normal(0,1, sample = True)

    ## SR shape =  N individuals---------------------------------------
    sr =  m.net.sender_receiver(
      focal_individual_predictors,
      target_individual_predictors,
      s_mu = 0.4, r_mu = -0.4
    )

    # Dyadic shape = N dyads--------------------------------------  
    dr = m.net.dyadic_effect(dyadic_predictors, d_sd=2.5) # Diadic effect intercept only 

    m.dist.bernoulli(logits = alpha + sr + dr, obs=network)

m.fit(model, num_samples = 500, num_warmup = 500, num_chains = 1, thinning = 1)

jax.local_device_count 32

  0%|          | 0/1000 [00:00<?, ?it/s]warmup:   0%|          | 1/1000 [00:01<20:22,  1.22s/it, 1 steps of size 2.34e+00. acc. prob=0.00]warmup:   1%|          | 10/1000 [00:01<01:38, 10.08it/s, 511 steps of size 1.25e-02. acc. prob=0.59]warmup:   2%|▏         | 16/1000 [00:01<01:02, 15.67it/s, 511 steps of size 2.08e-02. acc. prob=0.67]warmup:   2%|▏         | 24/1000 [00:01<00:38, 25.12it/s, 127 steps of size 5.01e-02. acc. prob=0.72]warmup:   3%|▎         | 34/1000 [00:01<00:25, 37.53it/s, 255 steps of size 3.06e-02. acc. prob=0.74]warmup:   4%|▍         | 42/1000 [00:01<00:22, 43.01it/s, 511 steps of size 1.22e-02. acc. prob=0.74]warmup:   5%|▌         | 50/1000 [00:01<00:18, 50.67it/s, 511 steps of size 9.93e-03. acc. prob=0.74]warmup:   6%|▌         | 58/1000 [00:02<00:17, 53.46it/s, 127 steps of size 3.73e-02. acc. prob=0.76]warmup:   6%|▋         | 65/1000 [00:02<00:16, 55.98it/s, 127 steps of size 3.43e-02. acc. prob=0.76]warmup:   7%|▋         | 72/1000 [00:02<00:16, 55.02it/s, 127 steps of size 1.31e-02. acc. prob=0.76]warmup:   8%|▊         | 80/1000 [00:02<00:15, 58.21it/s, 511 steps of size 1.76e-02. acc. prob=0.76]warmup:   9%|▉         | 88/1000 [00:02<00:14, 63.04it/s, 127 steps of size 4.08e-02. acc. prob=0.77]warmup:  10%|▉         | 97/1000 [00:02<00:13, 68.39it/s, 255 steps of size 3.27e-02. acc. prob=0.77]warmup:  11%|█         | 112/1000 [00:02<00:10, 88.69it/s, 127 steps of size 8.87e-02. acc. prob=0.77]warmup:  13%|█▎        | 131/1000 [00:02<00:07, 114.90it/s, 127 steps of size 5.36e-02. acc. prob=0.77]warmup:  14%|█▍        | 144/1000 [00:02<00:07, 117.95it/s, 255 steps of size 3.91e-02. acc. prob=0.77]warmup:  16%|█▌        | 157/1000 [00:03<00:06, 120.53it/s, 63 steps of size 1.21e-01. acc. prob=0.77] warmup:  18%|█▊        | 175/1000 [00:03<00:06, 137.48it/s, 31 steps of size 2.52e-01. acc. prob=0.78]warmup:  20%|██        | 204/1000 [00:03<00:04, 181.17it/s, 31 steps of size 1.24e-01. acc. prob=0.78]warmup:  23%|██▎       | 234/1000 [00:03<00:03, 215.59it/s, 31 steps of size 1.77e-01. acc. prob=0.78]warmup:  26%|██▌       | 260/1000 [00:03<00:03, 224.07it/s, 127 steps of size 4.67e-02. acc. prob=0.78]warmup:  28%|██▊       | 283/1000 [00:03<00:04, 168.74it/s, 63 steps of size 1.06e-01. acc. prob=0.78] warmup:  30%|███       | 303/1000 [00:03<00:04, 154.79it/s, 63 steps of size 2.25e-02. acc. prob=0.78]warmup:  32%|███▏      | 321/1000 [00:03<00:04, 146.99it/s, 127 steps of size 4.21e-02. acc. prob=0.78]warmup:  34%|███▎      | 337/1000 [00:04<00:04, 140.45it/s, 127 steps of size 4.74e-02. acc. prob=0.78]warmup:  36%|███▌      | 358/1000 [00:04<00:04, 156.25it/s, 63 steps of size 6.74e-02. acc. prob=0.78] warmup:  38%|███▊      | 375/1000 [00:04<00:04, 151.74it/s, 127 steps of size 5.04e-02. acc. prob=0.78]warmup:  41%|████      | 408/1000 [00:04<00:03, 195.12it/s, 63 steps of size 1.30e-01. acc. prob=0.78] warmup:  43%|████▎     | 434/1000 [00:04<00:02, 211.78it/s, 31 steps of size 1.72e-01. acc. prob=0.79]warmup:  46%|████▌     | 460/1000 [00:04<00:02, 223.16it/s, 63 steps of size 1.10e-01. acc. prob=0.78]warmup:  48%|████▊     | 484/1000 [00:04<00:02, 205.52it/s, 63 steps of size 1.23e-01. acc. prob=0.78]sample:  51%|█████     | 506/1000 [00:04<00:02, 195.53it/s, 63 steps of size 7.62e-02. acc. prob=0.90]sample:  53%|█████▎    | 527/1000 [00:04<00:02, 197.98it/s, 63 steps of size 7.62e-02. acc. prob=0.90]sample:  55%|█████▍    | 548/1000 [00:05<00:02, 195.92it/s, 63 steps of size 7.62e-02. acc. prob=0.90]sample:  57%|█████▋    | 568/1000 [00:05<00:02, 196.75it/s, 63 steps of size 7.62e-02. acc. prob=0.92]sample:  59%|█████▉    | 588/1000 [00:05<00:02, 195.38it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  61%|██████    | 610/1000 [00:05<00:01, 200.49it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  63%|██████▎   | 631/1000 [00:05<00:01, 203.06it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  65%|██████▌   | 652/1000 [00:05<00:01, 202.11it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  67%|██████▋   | 673/1000 [00:05<00:01, 197.94it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  69%|██████▉   | 693/1000 [00:05<00:01, 197.86it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  71%|███████▏  | 714/1000 [00:05<00:01, 199.28it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  73%|███████▎  | 734/1000 [00:06<00:01, 198.33it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  75%|███████▌  | 754/1000 [00:06<00:01, 190.33it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  77%|███████▋  | 774/1000 [00:06<00:01, 186.39it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  79%|███████▉  | 793/1000 [00:06<00:01, 185.75it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  81%|████████▏ | 814/1000 [00:06<00:00, 191.72it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  84%|████████▎ | 835/1000 [00:06<00:00, 193.41it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  86%|████████▌ | 856/1000 [00:06<00:00, 197.46it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  88%|████████▊ | 878/1000 [00:06<00:00, 202.08it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  90%|████████▉ | 899/1000 [00:06<00:00, 201.42it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  92%|█████████▏| 920/1000 [00:06<00:00, 202.50it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  94%|█████████▍| 943/1000 [00:07<00:00, 207.50it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  97%|█████████▋| 966/1000 [00:07<00:00, 211.53it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample:  99%|█████████▉| 988/1000 [00:07<00:00, 213.86it/s, 63 steps of size 7.62e-02. acc. prob=0.93]sample: 100%|██████████| 1000/1000 [00:07<00:00, 136.00it/s, 63 steps of size 7.62e-02. acc. prob=0.93]

R

library(BayesianInference)
# Setup platform------------------------------------------------
m=importBI(platform='cpu')

# Import data ------------------------------------------------
load(paste(system.file(package = "BayesianInference"),'/data/STRAND sim sr only.Rdata', sep = ''))

m$data_on_model = list()
m$data_on_model$N_id = length(ids)
m$data_on_model$network = m$net$mat_to_edgl(model_dat$outcomes[,,1])
m$data_on_model$N_dyads = m$net$mat_to_edgl(model_dat$outcomes[,,1])$shape[[1]]
m$data_on_model$focal_individual_predictors = jnp$array(model_dat$individual_predictors)
m$data_on_model$target_individual_predictors =jnp$array(model_dat$individual_predictors)

# Define model ------------------------------------------------
model <- function( N_id, N_dyads, network, focal_individual_predictors, target_individual_predictors){

  
  ## Block ---------------------------------------
  B = bi.dist.normal(0, 2.5, shape=c(1), name = 'block')
  
  #SR ---------------------------------------
  sr =  m$net$sender_receiver(focal_individual_predictors,target_individual_predictors)
  
  ### Dyadic--------------------------------------  
  #dr, dr_raw, dr_sigma, dr_L = m.net.dyadic_random_effects(idx.shape[0], cholesky_density = 2)# shape = n dyads
  dr = m$net$dyadic_effect(shape = c(N_dyads))
  
  ## SR ---------------------------------------                                                      
  m$dist$poisson(jnp$exp(B + sr + dr), obs=network)  
}

# Run MCMC ------------------------------------------------
m$fit(model) # Optimize model parameters through MCMC sampling

# Summary ------------------------------------------------
summary =m$summary()

Julia

# Setup device------------------------------------------------
using BayesianInference

# Setup device------------------------------------------------
m = importBI(platform="cpu")

# Simulate data ------------------------------------------------
N = 50
individual_predictor = m.dist.normal(0,1, shape = (N,1), sample = true)

kinship = m.dist.bernoulli(0.3, shape = (N,N), sample = true)
kinship = kinship.at[jnp.diag_indices(N)].set(0)

function sim_network(kinship, individual_predictor)
  # Intercept
  alpha = m.dist.normal(0,1, sample = true)

  # SR
  sr = m.net.sender_receiver(individual_predictor, individual_predictor, s_mu = 0.4, r_mu = -0.4, sample = true)

  # D
  DR = m.net.dyadic_effect(kinship, d_sd=2.5, sample = true)

  return m.dist.bernoulli(logits = alpha + sr + DR, sample = true)

end

network = sim_network(m.net.mat_to_edgl(kinship), individual_predictor)


# Predictive model ------------------------------------------------

m.data_on_model = pydict(
    network = network, 
    dyadic_predictors = m.net.mat_to_edgl(kinship),
    focal_individual_predictors = individual_predictor,
    target_individual_predictors = individual_predictor
)


@BI function model(network, dyadic_predictors, focal_individual_predictors, target_individual_predictors)
    N_id = network.shape[0]

    # Block ---------------------------------------
    alpha = m.dist.normal(0,1, name = "alpha")

    ## SR shape =  N individuals---------------------------------------
    sr =  m.net.sender_receiver(
      focal_individual_predictors,
      target_individual_predictors,
      s_mu = 0.4, r_mu = -0.4
    )

    # Dyadic shape = N dyads--------------------------------------  
    dr = m.net.dyadic_effect(dyadic_predictors, d_sd=2.5) # Diadic effect intercept only 

    m.dist.bernoulli(logits = alpha + sr + dr, obs=network)
end

m.fit(model, num_samples = 500, num_warmup = 500, num_chains = 1, thinning = 1)

Caution

Event if you don’t have dyadic effect, or block model effect, they need to be define to create intercepts (means) for those effects

Mathematical Details

Main Formula

The simple model that can be built to model link weights between nodes i and j can be defined using a Poisson distribution:

G_{ij} \sim \text{Poisson}(Y_{ij})

log(Y_{ij}) = \alpha + \lambda_i + \pi_j + \delta_{ij} + \beta_1 X_i + \beta_2 X_j + \beta_3 Q_{ij}

where:

• Y_{ij} is the weight of the link between i and j.

• \lambda_i is the sender random effect 🛈.

• \pi_j is the receiver random effect 🛈.

• \delta_{ij} is the dyadic random effect 🛈.

• \beta_1 is the effect of an individuals i level feature on the emission of a link (i.e., out-strength).

• \beta_2 is the effect of an individuals j level feature on the receiving a link (i.e., in-strength).

• \beta_3 is the effect of an dyadic characteristic between i and j on the likelihood of a tie.

Defining formula sub-equations and prior distributions

The sender and receiver random effects are similar to those described in chapter 13: Varying intercepts, but they are defined here using a joint prior so as to estimate the correlation within individuals to emit and receive a link:

\left(\begin{array}{cc} \lambda_i \\ \pi_i \end{array}\right) = \begin{array}{cc} \left(\begin{array}{cc} \sigma_\lambda \\ \sigma_\pi \end{array}\right) \circ \left(\begin{array}{cc} L \left(\begin{array}{cc} \hat{\lambda}_i \\ \hat{\pi}_i \end{array}\right) \end{array}\right) \end{array}

\sigma_\lambda \sim \text{Exponential}(1)

\sigma_\pi \sim \text{Exponential}(1)

L \sim \text{LKJ}(2)

\hat{\lambda}_i \sim \text{Normal}(0,1)

\hat{\pi}_i \sim \text{Normal}(0,1)

Similarly, for each dyad we can define a joint prior to estimate correlation between i–j links and j–i links:

\left(\begin{array}{cc} \delta_{ij} \\ \delta_{ji} \end{array}\right) = \begin{array}{cc} \left(\begin{array}{cc} \sigma_\delta \\ \sigma_\delta \end{array}\right) \circ \left(\begin{array}{cc} L_\delta \left(\begin{array}{cc} \hat{\delta}_{ij} \\ \hat{\delta}_{ji} \end{array}\right) \end{array}\right) \end{array}

\sigma_\delta \sim \text{Exponential}(1)

L_\delta \sim \text{LKJ}(2)

\hat{\delta}_{ij} \sim \text{Normal}(0,1)

Note(s)

Note

• Note that any additional covariates can be summed with a regression coefficient to \lambda_i, \pi_j and \delta_{ij}. Of course, for \lambda_i and \pi_j, as they represent nodal effects, these covariates need to be nodal characteristics (e.g., sex, age), whereas for \delta_{ij}, as it represents dyadic effects, these covariates need to be dyadic characteristics (e.g., genetic distances). Considering the previous example, given a vector of nodal characteristics, individual_predictors, and a matrix of dyadic characteristics, kinship, we can incorporate these covariates into the sender-receiver and dyadic effects, respectively.

• We can apply multiple variables as in chapter 2: Multiple Continuous Variables.

• We can apply interaction terms as in chapter 3: Interaction Between Continuous Variables.

• Network links can be modeled using Bernoulli (for proportions), Binomial (for unweighted network), Poisson or zero-inflated Poisson distributions (for count). In BI, you just need to set the correct likelihood distributions. For example, if you want to model the number of interactions between nodes, you can use the Poisson distribution. If you want to model the existence or absence of a link, you can use the Bernoulli distribution.

• If the network is undirected, then accounting for the correlation between the propensity to emit and receive links is not necessary, and the terms \lambda_i, \pi_j, and \delta_{ij} are no longer required. (Is it correct?)

• To account for exposure on a poisson model treat exposure as a nodal characteristic with its own parameter effect (i.e., regression coefficient). Their is several function that will help you to convert vectors or matrices in edge list format to have compatible data structure for the model (see API reference for bi.net.vec_to_edgl and bi.net.mat_to_edgl). F

• In the following chapters, we will see how to incorporate additional network effects into the model to account for network structural properties (e.g., clusters, assortativity, triadic closure, etc.).

Reference(s)

Ross, Cody T, Richard McElreath, and Daniel Redhead. 2024. “Modelling Animal Network Data in r Using STRAND.” Journal of Animal Ecology 93 (3): 254–66.